######## Replication file for "Talking while Fighting: Understanding the Role of Wartime Negotiation" ########

library(survival)
library(ggplot2)
library(stargazer)
library(lme4)
library(lmtest)
library(strucchange)
library(plyr)



setwd("./../figures")

##### IMPORT DATA #####

### War-day data
d = read.csv("./../replication/talkFightData.csv", stringsAsFactors = F)
d$date = as.Date(d$date)
d = d[with(d, order(cowNum, dayNum)),]

### Negotiation period data 
negPeriods = read.csv("./../replication/talkFightNegPeriods.csv", stringsAsFactors = F)

### Sources used for data
sources = read.csv("./../replication/talkFightSources.csv", stringsAsFactors = F)


##### TABLE 1: Information on sources #####

# Number of separate citations (410)
nrow(sources)

# Number of unique sources (355)
length(unique(sources$source)) 

# Table 1: Types of sources across all data
table(sources$sourceType)

# Table 1: Types of sources, disaggregated by pre-1945 and post-1945 wars 
table(sources$sourceType, sources$post45)


##### FIGURE 1: Four negotiation plots #####

### Figure 1a: Korean War
kw = d[d$warName=="Korean",]
qplot(kw$date, kw$negot, geom="line") + geom_line(size=2) +
  ylab("Negotiation") + xlab("Date") + 
  scale_y_continuous(breaks=c(0,1), labels=c("No", "Yes"), limits=c(-0.2, 1.2)) + theme_bw() +
  scale_x_date(date_breaks="6 months", date_labels="%b %Y")

ggsave("korean.pdf", height=2, width=6)
rm(kw)

### Figure 1b: Roman Republic 
rr = d[d$warName=="Roman Republic",]
qplot(rr$date, rr$negot, geom="line") + geom_line(size=2) +
  ylab("Negotiation") + xlab("Date") +  
  scale_y_continuous(breaks=c(0,1), labels=c("No", "Yes"), limits=c(-0.2, 1.2)) + theme_bw() + 
  scale_x_date(date_breaks="1 months", date_labels="%b %Y")

ggsave("roman.pdf", height=2, width=6)
rm(rr)

### Figure 1c: Fourth Central American 
fc = d[d$warName=="Fourth Central American",]
qplot(fc$date, fc$negot, geom="line") + geom_line(size=2) +
  ylab("Negotiation") + xlab("Date") + 
  scale_y_continuous(breaks=c(0,1), labels=c("No", "Yes"), limits=c(-0.2, 1.2)) + theme_bw() +
  scale_x_date(date_breaks="2 weeks", date_labels="%b %d %Y")

ggsave("fourth.pdf", height=2, width=6)
rm(fc)

### Figure 1d: Ugandan-Tanzanian 
ut = d[d$warName=="Ugandian-Tanzanian",]
qplot(ut$date, ut$negot, geom="line") + geom_line(size=2) +
  ylab("Negotiation") + xlab("Date") + 
  scale_y_continuous(breaks=c(0,1), labels=c("No", "Yes"), limits=c(-0.2, 1.2)) + theme_bw() +
  scale_x_date(date_breaks="1 months", date_labels="%b %Y")

ggsave("uganda.pdf", height=2, width=6)
rm(ut)


##### FIGURE 2: Descriptive statistics #####

### Calculate statistics
fr = list(); wars = unique(d$warName)

for (i in 1:length(wars)) {
  propNeg = propFirst = post45 = NA
  
  # Get a war 
  oneWar = d[d$warName==wars[i],]
  
  # Calculate proportion of war with negotiation
  propNeg = sum(oneWar$negot)/nrow(oneWar) 
  propNeg = ifelse(propNeg==0, -0.05, propNeg)
  
  # Calculate proportion of war elapsed before first negotiation
  propFirst = min(which(oneWar$negot==1))/nrow(oneWar)
  propFirst = ifelse(propFirst==Inf, 1.05, propFirst)
  post45 = oneWar$post45[1]

  # Collect data
  first = data.frame(warName=oneWar$warName[1], post45=post45, propNeg=propNeg, propFirst=propFirst)
  fr[[i]] = first
}

firsts = do.call("rbind", fr)
rm(fr, first, propNeg, propFirst, oneWar, post45)


### Figure 2a: Proportion of war with ongoing negotiations 
ggplot(firsts, aes(x=propNeg)) + geom_histogram(aes(fill=propNeg<0), binwidth=0.05) + theme_bw() + ylab("Count") + 
  xlab("Proportion of War with Ongoing Negotiations") + ylim(0,20) +
  scale_fill_manual(values=c("gray30", "gray70"), guide=guide_legend(reverse=T)) + 
  guides(fill=FALSE) + scale_x_continuous(breaks=c(0, 0.25, 0.5, 0.75, 1)) 
ggsave("propWithNeg.pdf", height=3, width=4)


### Figure 2b: Proportion of war elapsed before first negotiation 
ggplot(firsts, aes(x=propFirst)) + geom_histogram(aes(fill=propFirst>1), binwidth=0.075) + theme_bw() + ylab("Count") + 
  xlab("Proportion of War Elapsed before First Negotiation") +
  scale_fill_manual(values=c("gray30", "gray70"), guide=guide_legend(reverse=T)) + guides(fill=FALSE) + scale_x_continuous(breaks=c(0, 0.25, 0.5, 0.75, 1)) +
  ylim(0,20) 
ggsave("propFirstNeg.pdf", height=3, width=4)


### Figure 2c: Total number of negotiation days
# Calculate number of negotiation periods and total number of negotiation days
numNegs = negDays = NA
for (i in 1:length(wars)) {
  oneWar = negPeriods[which(negPeriods$warName==wars[i]),]
  numNegs[i] = nrow(oneWar)
  negDays[i] = sum(oneWar$negLength)
}

qplot(numNegs) + geom_histogram(binwidth=1, aes(fill=numNegs==0)) +
  xlab("Number of Negotiation Periods") + ylab("Count") + scale_x_continuous(breaks=c(0,2,4,6,8,10)) +
  scale_fill_manual(values=c("gray30", "gray70")) + theme_bw() + guides(fill=FALSE) + ylim(0,30)
ggsave("numNegHist.pdf", height=3, width=4)


### Figure 2d: Lengths of negotiations 
qplot(negPeriods$negLength[which(negPeriods$negLength < 150)]) + geom_histogram(binwidth=10) +
  scale_x_continuous(breaks=seq(0,150,25)) +
  scale_y_continuous(breaks=seq(0, 80, 20), limits = c(0,80)) + 
  xlab("Negotiation Period Length (Days)") + ylab("Count") +
  guides(fill=FALSE) + theme_bw() 
ggsave("negLenHist.pdf", height=3, width=4)


##### TABLE 2: Distribution of the negotiation variable ####

table(d$negot)
prop.table(table(d$negot))

table(d$post45, d$negot)
prop.table(table(d$post45, d$negot), margin=1)


##### IN-TEXT: Structural break test #####

yr = 1820:2005 
nw = nwd = ng = NA

for (i in 1:length(yr)) {
  oy = d[format(d$date, "%Y")==yr[i],]
  nw[i] = length(unique(oy$warName))   # Number of wars
  nwd[i] = nrow(oy)                    # Number of war-days
  ng[i] = sum(oy$negot)                # Number of negotiation days
}

# Create time series
negs = ts(ng, start=1820)

# Identify break points
breaktest = breakpoints(negs ~ nw + nwd)
summary(breaktest) # Best model in BIC has two cutpoints at 1945 and 1972

rm(oy, yr, nw, nwd, ng)


##### IN-TEXT: Logistic regressions #####

y = 1824:2002
bics = est = se = NA

for (i in 1:length(y)) {
  ntest = glm(negot ~ I(format(date, "%Y") > y[i]), data=d, family=binomial)
  bics[i] = BIC(ntest)
  est[i] = summary(ntest)$coefficients[2,1]
  se[i] = summary(ntest)$coefficients[2,2]
}
bicout = data.frame(year=y, bic=bics, est, se)

# Identify year(s) with lowest BIC
bicout$year[which(bicout$bic==min(bicout$bic))]


rm(y, bics, est, se)



##### FIGURE 3: Boxplots of negotiation patterns #####

### Figure 3a: Proportion of war with ongoing negotiations 
qplot(factor(firsts$post45[which(firsts$propNeg>0)], levels=c(1,0)), firsts$propNeg[which(firsts$propNeg>0)], geom="boxplot") + theme_bw() + coord_flip() +
  scale_x_discrete(name="Era", labels=c("Post-1945", "Pre-1945")) +
  scale_y_continuous(name="Proportion of War with Ongoing Negotiations", breaks=c(0, 0.25, 0.5, 0.75, 1), limits=c(0,1)) 
ggsave("propNegBox.pdf", height=2, width=5)

### Figure 3b: Proportion of war elapsed before first negotiation 
qplot(factor(firsts$post45[which(firsts$propFirst <= 1)], levels=c(1,0)), firsts$propFirst[which(firsts$propFirst <= 1)], geom="boxplot") + theme_bw() + coord_flip() +
  scale_x_discrete(name="Era", labels=c("Post-1945", "Pre-1945")) +
  scale_y_continuous(name="Proportion of War Elapsed before First Negotiation", breaks=c(0, 0.25, 0.5, 0.75, 1), limits=c(0,1)) 
ggsave("propFirstBox.pdf", height=2, width=5)



##### FIGURE 4: Negotiation trajectories #####

### Create necessary trajectory data
pre = d[which(d$post45==0),]
post = d[which(d$post45==1),]
warsPre = unique(pre$warName)
warsPost = unique(post$warName)

rseq = seq(0, 1, 0.05)
meansPre = meansPost = list()

# Pre-1945 trajectories
for (i in 1:length(warsPre)) {
  oneWar = d[d$warName==warsPre[i],]
  mres = NA
  for (j in 1:length(rseq)) {
    part = oneWar[oneWar$propDoneR==rseq[j],]
    mres[j] = mean(part$negot, na.rm=T)
  }
  meansPre[[i]] = data.frame(propDone=rseq, negot=mres)
}
mrPre = do.call('rbind', meansPre)
mrPre$negot = na.approx(mrPre$negot, na.rm=F); 
mrPre = na.omit(mrPre)
smPre = ksmooth(mrPre$propDone, mrPre$negot, bandwidth=0.2)

trajPre = data.frame(propDone=smPre$x, negot=smPre$y)
trajPre$propDone = round_any(trajPre$propDone, 0.05, floor)
trajPre = unique(trajPre)
trajPre = trajPre[!duplicated(trajPre$propDone),]
trajPre$era = "Pre-1945"

# Post-1945 trajectories
for (i in 1:length(warsPost)) {
  oneWar = d[d$warName==warsPost[i],]
  mres = NA
  for (j in 1:length(rseq)) {
    part = oneWar[oneWar$propDoneR==rseq[j],]
    mres[j] = mean(part$negot, na.rm=T)
  }
  meansPost[[i]] = data.frame(propDone=rseq, negot=mres)
}
mrPost = do.call('rbind', meansPost)
mrPost$negot = na.approx(mrPost$negot, na.rm=F); 
mrPost = na.omit(mrPost)
smPost = ksmooth(mrPost$propDone, mrPost$negot, bandwidth=0.2)

trajPost = data.frame(propDone=smPost$x, negot=smPost$y)
trajPost$propDone = round_any(trajPost$propDone, 0.05, floor)
trajPost = unique(trajPost)
trajPost = trajPost[!duplicated(trajPost$propDone),]
trajPost$era = "Post-1945"

rm(part, mres, oneWar, mrPre, mrPost, smPre, smPost, rseq)

# Put pre-1945 and post-1945 trajectory datasets together
traj = rbind(trajPre, trajPost); rm(trajPre, trajPost)
traj$era = factor(traj$era, levels=c("Pre-1945", "Post-1945"))

### Figure 4: Prevalence of negotiations over the course of war
ggplot(traj, aes(x=propDone, y=negot, group=era)) + geom_line(aes(color=era), lineend="round", size=1.5) + theme_bw() + 
  xlab("Proportion of War Elapsed") + ylab("Pr(Negotiation)") +
  scale_color_manual("", values=c("black", "gray60"), labels=c("Pre-1945", "Post-1945")) + ylim(0,0.45) 
ggsave("smoothNegPP.pdf", height=3, width=5.5)


##### TABLE 3: Cox proportional hazard models #####
##### FIGURE 5: Marginal effects plot #####

### Setting up survival analysis 
surv = Surv(as.numeric(d$dayNum), as.numeric(d$dayNumNext), as.numeric(d$terminate))

# Model 1: Negotiations, bivariate
pFit1a = coxph(surv ~ negot + cluster(warName), data=d)
summary(pFit1a)
cox.zph(pFit1a)

# Model 2: Negotiations, controls
pFit1b = coxph(surv ~ negot + issues + contig + cincRatio + demo + majNuk + 
                 issues:logDayNum + cluster(warName), data=d)
summary(pFit1b)
cox.zph(pFit1b)

# Model 3: Negotiations x Post-1945, bivariate
pFit1c = coxph(surv ~ post45*negot + 
                 post45:logDayNum + post45:negot:logDayNum + cluster(warName), data=d)
summary(pFit1c)
cox.zph(pFit1c)

# Model 4: Negotiations x Post-1945, controls
pFit1d = coxph(surv ~ post45*negot + issues + contig + cincRatio + demo + majNuk +
                 post45:logDayNum + cluster(warName), data=d)
summary(pFit1d)
cox.zph(pFit1d)


### Table 3 (Remove this before submitting)
stargazer(pFit1a, pFit1b, pFit1c, pFit1d, 
          covariate.labels = c("Post-1945", "Negotiation", "Issue salience", "Contiguity",
                               "CINC ratio", "Democracy", "Major/Nuclear", "Remove",
                               "Post-1945 $\times$ Negotiation", "Remove", "Remove"),
          column.sep.width = "-5pt", no.space=T, align=T, df=F)



### Calculate marginal rates (see Brambor et al. 2006)
cv = vcov(pFit1d)
v45 = cv["post45", "post45"]
vNeg = cv["negot", "negot"]
vRF = cv["post45:negot", "post45:negot"]
covFI = cv["post45:negot","negot"]
negCoef = coef(summary(pFit1d))["negot", "coef"]
negCoefSE = coef(summary(pFit1d))["negot", "se(coef)"]
p45Coef = coef(summary(pFit1d))["negot", "coef"]
p45CoefSE = coef(summary(pFit1d))["post45", "se(coef)"]
rfCoef = coef(summary(pFit1d))["post45:negot", "coef"]

se0 = sqrt(vNeg + 0^2*vRF + 2*0*covFI)
se1 = sqrt(vNeg + 1^2*vRF + 2*1*covFI)
ses = data.frame(x=c("Pre-1945", "Post-1945"), y=c(negCoef, negCoef+rfCoef), se=c(se0, se1))
ses$x = factor(ses$x)


### Figure 5: Marginal effect of negotiations on the termination of conflict
ggplot(ses, aes(x=x, y=y)) + geom_point(size=3, color="gray40") + geom_hline(yintercept=0, lty=2) + 
  geom_errorbar(aes(ymin=y-1.96*se, ymax=y+1.96*se), width=0.14, color="gray40", size=1) + 
  theme_bw() + coord_flip() + xlab("Era") + ylab("Marginal Effect of Negotiation on War Termination") 
ggsave("margEffectPP.pdf", height=2.25, width=5.5)

rm(cv, v45, vNeg, vRF, covFI, negCoef, negCoefSE, p45Coef, p45CoefSE, rfCoef, se0, se1)


##### TABLE 4: Negotiations and third-party pressure #####

### Table 4 itself
table(negPeriods$post45, negPeriods$pressure) 

### Chi-squared test
chisq.test(negPeriods$post45, negPeriods$pressure) 


##### FIGURE 6: Coefficient plots of potential mechanisms #####

### Identify potential mechanisms 
controls = c("none", "unsc", "logAllies", "nuk", "lawsOfWar", "ratif", "territ", "conquest")

est45s = se45s = estcons = secons = NA

for (i in 1:length(controls)) {
  
  if (i==1) {  # No mechanism, just negotiations
    form = as.formula("negot ~ negot1 + post45 + dayNum")
    ntest = glm(form, data=d, family=binomial, control = list(maxit = 100))
    est45 = coef(summary(ntest))["post45","Estimate"]
    se45 = coef(summary(ntest))["post45","Std. Error"]
    estcon = 0
    secon = 0
  } else {    # Add a mechanism
    form = as.formula(paste("negot ~ negot1 + post45 + dayNum +", controls[i], sep=" "))
    ntest = glm(form, data=d, family=binomial, control = list(maxit = 100))
    est45 = coef(summary(ntest))["post45","Estimate"]
    se45 = coef(summary(ntest))["post45","Std. Error"]
    estcon = coef(summary(ntest))[controls[i],"Estimate"]
    secon = coef(summary(ntest))[controls[i],"Std. Error"]
  }
  est45s[i] = est45
  se45s[i] = se45
  estcons[i] = estcon
  secons[i] = secon
}

mechOutput = data.frame(control=c(controls, controls), coef = c(est45s, estcons), se = c(se45s, secons),
                    Variable=c(rep("Post-1945", length(controls)), rep("Variable", length(controls))))
mechOutput$control = factor(mechOutput$control, levels = controls)
mechOutput$none = factor(ifelse(mechOutput$control=="none" & mechOutput$Variable=="Variable", 1, 0))


### Figure 6: Estimated coefficients for each potential mechanism
ggplot(mechOutput, aes(control, coef, colour=Variable, alpha=none)) +
  geom_hline(yintercept=0, linetype=2) + 
  geom_pointrange(aes(ymin = coef-1.96*se, ymax = coef+1.96*se), position=position_dodge(width=0.4)) +
  scale_x_discrete("Variable", labels=c("None", "UNSC", "Alliances", "Nuclear", "Laws of War", "Ratification", "Territorial", "Conquest")) +
  scale_alpha_manual(guide=FALSE, breaks=c(1,0), values=c(1,0)) +
  theme_bw() + scale_y_continuous("Coefficient") + scale_color_manual("", values=c("gray20", "gray60"))
ggsave("coefPlot45.pdf", height=3, width=8)

rm(est45, se45, est45s, se45s, estcon, secon, estcons, secons, form)



######## APPENDIX MATERIALS #########

##### APPENDIX FIGURE 1: Negotiation variable for three wars #####
# (See code above for "FIGURE 1: Four negotiation plots.")


##### APPENDIX TABLE 1: Negotiations by war #####
cbind(wars, numNegs, negDays)


##### APPENDIX TABLE 2: Descriptive statistics on wars and negotiations #####

# War lengths
summary(d$dayNum[which(d$terminate==1)])
summary(pre$dayNum[which(pre$terminate==1)])
summary(post$dayNum[which(post$terminate==1)])

# Number of negotiation periods per war
summary(numNegs)
summary(numNegs[1:57])
summary(numNegs[58:92])

# Negotiation period lengths
summary(negPeriods$negLength)
summary(negPeriods$negLength[which(negPeriods$post45==0)])
summary(negPeriods$negLength[which(negPeriods$post45==1)])

# Proportion of war with negotiation
summary(firsts$propNeg[which(firsts$propNeg >= 0)])
summary(firsts$propNeg[which(firsts$propNeg >= 0 & firsts$post45 == 0)])
summary(firsts$propNeg[which(firsts$propNeg >= 0 & firsts$post45 == 1)])

# Proportion of war elapsed before first negotiation 
summary(firsts$propFirst[which(firsts$propFirst <= 1)])
summary(firsts$propFirst[which(firsts$propFirst <= 1 & firsts$post45 == 0)])
summary(firsts$propFirst[which(firsts$propFirst <= 1 & firsts$post45 == 1)])


##### APPENDIX TABLE 3: Descriptive statistics for continuous variables #####
summary(d$issues); summary(pre$issues); summary(post$issues)
summary(d$cincRatio); summary(pre$cincRatio); summary(post$cincRatio)


##### APPENDIX TABLE 4: Descriptive statistics for binary variables #####
table(d$negot); table(pre$negot); table(post$negot)
table(d$contig); table(pre$contig); table(post$contig)
table(d$demo); table(pre$demo); table(post$demo)
table(d$majNuk); table(pre$majNuk); table(post$majNuk)




##### APPENDIX TABLE 5: Structural break test results #####
# (First run code above for "IN-TEXT: Structural break test.")

summary(breaktest) # Best model in BIC has two cutpoints at 1945 and 1972


##### APPENDIX TABLE 6: 95% confidence intervals for structural break tests #####

confint(breaktest, 1)
confint(breaktest, 2)
confint(breaktest, 3)
confint(breaktest, 4)
confint(breaktest, 5) 


##### APPENDIX FIGURE 2: BICs for logistic regressions #####
# (First run code above for "IN-TEXT: Logistic regressions.")

ggplot(bicout, aes(year, bic)) + theme_bw() + xlab("Cutoff Year") + ylab("BIC") + 
  geom_vline(xintercept=1945, color="gray50", size=0.7, linetype=2) +
  geom_line(color="gray20", size=1) 
ggsave("negBIC.pdf", height=3, width=5)


##### APPENDIX TABLE 7: Regressions regarding post-1945 impact on negotiation frequency #####

# Logistic with lagged DV, bivariate
ntest1a = glm(negot ~ negot1 + post45, data=d, family=binomial, control = list(maxit = 100))
summary(ntest1a)
ntest1ar = coeftest(ntest1a, vcov=vcovHC(ntest1a, type="HC3", cluster="warName", adjust=T))
ntest1ar

# Logistic with lagged DV, controls
ntest1b = glm(negot ~ negot1 + post45 + issues + contig + cincRatio + demo + majNuk + logDayNum,
              data=d, family=binomial)
summary(ntest1b)
ntest1br = coeftest(ntest1b, vcov=vcovHC(ntest1b, type="HC3", cluster="warName", adjust=T))
ntest1br

# Mixed effects OLS, bivariate
mtest1a = lmer(negot ~ post45 + (1|warName), data=d)
summary(mtest1a)

# Mixed effects OLS, controls
mtest1b = lmer(negot ~ post45 + issues + contig + cincRatio + demo + majNuk + logDayNum + (1|warName), data=d)
summary(mtest1b)


### Table 7: Results

# Create base table
stargazer(ntest1a, ntest1b, mtest1a, mtest1b, 
          covariate.labels = c("Lagged DV", "Post-1945", "Issue salience", "Contiguity",
                               "CINC ratio", "Democracy", "Major/Nuclear", "Time trend"),
          column.sep.width = "-5pt", no.space=T, align=T, df=F)

# Results with clustered SEs for logistic regressions (actual results in Table 7)
stargazer(ntest1ar, ntest1br, mtest1a, mtest1b,         
          covariate.labels = c("Lagged DV", "Post-1945", "Issue salience", "Contiguity",
                               "CINC ratio", "Democracy", "Major/Nuclear", "Time trend"),
          column.sep.width = "-5pt", no.space=T, align=T, df=F)
